
require(stats)
require(graphics)
#library("ggplot2")
require(ggplot2)
#library("grDevices")
require(grDevices)
#library("gridExtra")
require(gridExtra)
#library("grid")
require(grid)
#install.packages("Hmisc")
require(Hmisc)

# Please make sure to change the path to where the appropriate dataset is located on your computer

# LIFECYCLE
# FR
fr.k1.data <- read.csv("PNAS_FULL_RESTRICT_K1_LC.csv")
# ST
st.k1.data <- read.csv("PNAS_IND_ST_RESTRICT_K1_LC.csv")
# FAM
fam.k1.data <- read.csv("PNAS_IND_FAM_RESTRICT_K1_LC.csv")
#LS
ls.k1.data <- read.csv("PNAS_IND_LS_RESTRICT_K1_LC.csv")
#HS
hs.k1.data <- read.csv("PNAS_IND_HS_RESTRICT_K1_LC.csv")
#O
open.k1.data <- read.csv("PNAS_OPEN_K1_LC.csv")
#NO
no.k1.data <- read.csv("PNAS_NOCHAN_OPEN_K1_LC.csv")

# ORIGINAL 
# FULL RESTRICT
fr.k1.data1 <- read.csv("PNAS_VARILL_FR_K1.csv")
# STUDENT
st.k1.data1 <- read.csv("PNAS_VARILL_ST_K1.csv")
# FAMILY
fam.k1.data1 <- read.csv("PNAS_VARILL_FAM_K1.csv")
#LOW SKILLED
ls.k1.data1 <- read.csv("PNAS_VARILL_LS_K1.csv")
#HIGH SKILLED
hs.k1.data1 <- read.csv("PNAS_VARILL_HS_K1.csv")
#BASELINE
open.k1.data1 <- read.csv("PNAS_VARILL_O_K1.csv")
#FREE MOVEMENT
no.k1.data1 <- read.csv("PNAS_VARILL_NO_K1.csv")


sponsor.require <- 1
high.skilled.p <-  2
low.skilled.p <-  3
p.tourist.visa <- 4
p.fakedocs <-  5
no.econ.stud <- 6
p.all.illegal <- 7
k <- 8
my.seed <- 9

legaltourist.attempts <- c(10:29)

ability.legal.cols <- c(30:90)
ability.fakedocs.cols <- c(91:151)
ability.work.cols <- c(152:212)

migration.rate.cols <- c(213:273)
legal.cols <- c(274:334)
fi.cols <- c(335:395)
sl.cols <- c(396:456)
stud.mig <- c(457:517)
hs.mig <- c(518:578)
ls.mig <- c(579:639)
fam.mig <- c(640:700)
aspire.cols <- (701:761)

num.agents <-  1166 
aspire.enough <- 507
select.policy.level <-  0.3 #1-0.7

# df.list1 = list of data for original model

df.list1 <- list(no.k1.data1,fr.k1.data1,fam.k1.data1,ls.k1.data1,open.k1.data1,hs.k1.data1,st.k1.data1)

table(fr.k1.data1[,p.all.illegal])
table(st.k1.data1[,p.all.illegal])
table(fam.k1.data1[,p.all.illegal])
table(ls.k1.data1[,p.all.illegal])
table(hs.k1.data1[,p.all.illegal])
table(open.k1.data1[,p.all.illegal])
table(no.k1.data1[,p.all.illegal])

# Reduce "original" dataset to include only the data with the enforcement values (p.all.illegal) 
# shown in Figures 2-4. 

df.list1 <- lapply(df.list1, function(x) {
  x <- x[which(x[,p.all.illegal] == select.policy.level),]
  x} )

# df.list = list of data for lifecycle model. The lifecycle model (LC) uses another variable to control 
# enforcement: p.tourist.visa and p.fakedocs, not p.all.illegal. The  LC model does not have multiple 
#settings for enforcement (only that shown in Figs-2-4) and, therefore, lines just above do not apply.

df.list <- list(no.k1.data,fr.k1.data,fam.k1.data,ls.k1.data,open.k1.data,hs.k1.data,st.k1.data)

# Reduce each dataset to have the same number of rows for subtraction
df.list <- lapply(df.list, function(x) {
  x <- x[1:980,]
  x} )

df.list1 <- lapply(df.list1, function(x) {
  x <- x[1:980,]
  x} )

# Name variables
names.df.list <-  Cs(no.k1.data,fr.k1.data,fam.k1.data,ls.k1.data,open.k1.data,hs.k1.data,st.k1.data)

names.df.list1 <-  Cs(no.k1.data1,fr.k1.data1,fam.k1.data1,ls.k1.data1,open.k1.data1,hs.k1.data1,st.k1.data1)

names(df.list) <- names.df.list
names(df.list1) <- names.df.list1

time.point <- 20

# Measure variables at time.point

#Lifecycle
for(i in 1:length(df.list))
{
  df.list[[i]]$legal.raw <- (df.list[[i]][,legal.cols[time.point]])
  df.list[[i]]$legal <-  df.list[[i]]$legal.raw / num.agents
  df.list[[i]]$fill <- (df.list[[i]][,fi.cols[time.point]])
  df.list[[i]]$semi <- (df.list[[i]][,sl.cols[time.point]])
  df.list[[i]]$all.illegal.raw <- (df.list[[i]]$fill + df.list[[i]]$semi)
  df.list[[i]]$all.illegal <- df.list[[i]]$all.illegal.raw / num.agents
  df.list[[i]]$nm <- 1 - (df.list[[i]]$legal + df.list[[i]]$all.illegal) # nonmigrants
  df.list[[i]]$all.migs.raw <- df.list[[i]]$legal.raw + df.list[[i]]$all.illegal.raw
}

# Original
for(i in 1:length(df.list1))
{
  df.list1[[i]]$legal.raw <- (df.list1[[i]][,legal.cols[time.point]])
  df.list1[[i]]$legal <-  df.list1[[i]]$legal.raw / num.agents
  df.list1[[i]]$fill <- (df.list1[[i]][,fi.cols[time.point]])
  df.list1[[i]]$semi <- (df.list1[[i]][,sl.cols[time.point]])
  df.list1[[i]]$all.illegal.raw <- (df.list1[[i]]$fill + df.list1[[i]]$semi)
  df.list1[[i]]$all.illegal <- df.list1[[i]]$all.illegal.raw / num.agents
  df.list1[[i]]$nm <- 1 - (df.list1[[i]]$legal + df.list1[[i]]$all.illegal) # nonmigrants
  df.list1[[i]]$all.migs.raw <- df.list1[[i]]$legal.raw + df.list1[[i]]$all.illegal.raw
}


  
percentile <- function(x){
  quantile(x,c(0.0275, 0.975))
}

# SUBTRACT BASELINE

# Lifecycle

# Find index for baseline
b.i <- grep("open", names.df.list)

for (i in 1:7)
{
  df.list[[i]]$legal.adj <- df.list[[i]]$legal - df.list[[b.i]]$legal
  df.list[[i]]$nm.adj  <- df.list[[i]]$nm - df.list[[b.i]]$nm
  df.list[[i]]$illegal.adj  <- df.list[[i]]$all.illegal - df.list[[b.i]]$all.illegal
  df.list[[i]]$legal.migs.adj <- (df.list[[i]]$legal.raw - df.list[[b.i]]$legal.raw) / df.list[[i]]$all.migs.raw
  df.list[[i]]$illegal.migs.adj  <- (df.list[[i]]$all.illegal.raw - df.list[[b.i]]$all.illegal.raw) / df.list[[i]]$all.migs.raw
}

# Original

b1.i <- grep("open", names.df.list1)

for (i in 1:7)
{
  df.list1[[i]]$legal.adj <- df.list1[[i]]$legal - df.list1[[b1.i]]$legal
  df.list1[[i]]$nm.adj  <- df.list1[[i]]$nm - df.list1[[b1.i]]$nm
  df.list1[[i]]$illegal.adj  <- df.list1[[i]]$all.illegal - df.list1[[b1.i]]$all.illegal
  df.list1[[i]]$legal.migs.adj <- (df.list1[[i]]$legal.raw - df.list1[[b1.i]]$legal.raw) / df.list1[[i]]$all.migs.raw
  df.list1[[i]]$illegal.migs.adj  <- (df.list1[[i]]$all.illegal.raw - df.list1[[b1.i]]$all.illegal.raw) / df.list1[[i]]$all.migs.raw
}

# MAKE SMALLER VECTORS TO POPULATE WITH MEAN, LOWER BOUND AND UPPER BOUND

lc <-  vector("list", length = 3)
or <-  vector("list", length = 3)

names(lc) <- Cs(mean.df, lb, ub)
names(or) <- Cs(mean.df, lb, ub)

names.temp <- c("Free Movement","Closed","Family","Low-Skilled","Baseline","High-Skilled","Student")

for (i in 1: length(lc))
{
lc[[i]] <- matrix(data = NA, nrow = 7, ncol = 3)
colnames(lc[[i]]) <- Cs(leg,ill,nm)
rownames(lc[[i]]) <- names.temp 

or[[i]] <- matrix(data = NA, nrow = 7, ncol = 3)
colnames(or[[i]]) <- Cs(leg,ill,nm)
rownames(or[[i]]) <- names.temp 
}

# POPULATE 

#MEANS

for (i in 1:length(df.list))
{
  lc[[1]][i,1] <- mean(df.list[[i]]$legal.adj)
  lc[[1]][i,2] <- mean(df.list[[i]]$illegal.adj)
  lc[[1]][i,3] <- mean(df.list[[i]]$nm.adj)
  or[[1]][i,1] <- mean(df.list1[[i]]$legal.adj)
  or[[1]][i,2] <- mean(df.list1[[i]]$illegal.adj)
  or[[1]][i,3] <- mean(df.list1[[i]]$nm.adj)
}

# LB

for (i in 1:length(df.list))
{
  lc[[2]][i,1] <- percentile(df.list[[i]]$legal.adj)[1]
  lc[[2]][i,2] <- percentile(df.list[[i]]$illegal.adj)[1]
  lc[[2]][i,3] <- percentile(df.list[[i]]$nm.adj)[1] 
  or[[2]][i,1] <- percentile(df.list1[[i]]$legal.adj)[1]
  or[[2]][i,2] <- percentile(df.list1[[i]]$illegal.adj)[1]
  or[[2]][i,3] <- percentile(df.list1[[i]]$nm.adj)[1] 
}

# UB

for (i in 1:length(df.list))
{
  lc[[3]][i,1] <- percentile(df.list[[i]]$legal.adj)[2]
  lc[[3]][i,2] <- percentile(df.list[[i]]$illegal.adj)[2]
  lc[[3]][i,3] <- percentile(df.list[[i]]$nm.adj)[2] 
  or[[3]][i,1] <- percentile(df.list1[[i]]$legal.adj)[2]
  or[[3]][i,2] <- percentile(df.list1[[i]]$illegal.adj)[2]
  or[[3]][i,3] <- percentile(df.list1[[i]]$nm.adj)[2] 
}


# make percentage 

for (i in 1:length(lc))
{
  lc[[i]] <- lc[[i]]*100
  or[[i]] <- or[[i]]*100
}

#Remove Baseline (=0)

for (i in 1:length(lc))
{
  lc[[i]] <- lc[[i]][-5,]
  or[[i]] <- or[[i]][-5,]
}

plot.obj.names <- c("Free Movement","Closed","Family","Low-Skilled","High-Skilled","Student")

num.sets <- 1:nrow(lc$mean.df)
jitter<-.2

pdf("FigF1_coeff.pdf", width = 7, height = 7, onefile = FALSE)
lay <- layout(matrix(c(1,2,3)),heights = c(4,4,4), widths = c(5,5,5))
par(mai =c(0.6,2,0.3,2), cex.axis=1.5, xpd = TRUE)
plot(lc$mean.df[,1], num.sets, type="n", ylab ="",xlab ="", yaxt="n", xlim = c(-23,23), ylim = c(0.8,6.2))
segments(lc$lb[,1], num.sets, lc$ub[,1], num.sets, col = "purple3")
points(lc$mean.df[,1], num.sets, pch= 19, cex=1, col = "purple3")
segments(or$lb[,1], num.sets+jitter, or$ub[,1], num.sets+jitter, col = "black")
points(or$mean.df[,1], num.sets+jitter, pch= 19, cex=1, col = "black")
mtext("Legal",side=3,line=1,outer=F, font = 2)
axis(2, at=num.sets, labels=plot.obj.names, las = 2)
legend(25.5, 6, pch=c(19,19), legend=c("Life Cycle", "Original"), col = c("purple3", "black"),
       cex =1.2, border = "black", horiz = F)

par( mai =c(0.6,2,0.3,2), cex.axis=1.5, xpd = TRUE)
plot(lc$mean.df[,2], num.sets, type="n", ylab ="",xlab ="", yaxt="n", xlim = c(-22,22), ylim = c(0.8,6.2))
segments(lc$lb[,2], num.sets, lc$ub[,2], num.sets, col = "purple3")
points(lc$mean.df[,2], num.sets, pch= 19, cex=1, col = "purple3")
segments(or$lb[,2], num.sets+jitter, or$ub[,2], num.sets+jitter, col = "black")
points(or$mean.df[,2], num.sets+jitter, pch= 19, cex=1, col = "black")
mtext("Unauthorized",side=3,line=1,outer=F, font = 2)
axis(2, at=num.sets, labels=plot.obj.names, las = 2)

par( mai =c(0.6,2,0.3,2), cex.axis=1.5, xpd = TRUE)
plot(lc$mean.df[,3], num.sets, type="n", ylab ="", xlab ="", yaxt="n", xlim = c(-22,22), ylim = c(0.8,6.2))
segments(lc$lb[,3], num.sets, lc$ub[,3], num.sets, col = "purple3")
points(lc$mean.df[,3], num.sets, pch= 19, cex=1, col = "purple3")
segments(or$lb[,3], num.sets+jitter, or$ub[,3], num.sets+jitter, col = "black")
points(or$mean.df[,3], num.sets+jitter, pch= 19, cex=1, col = "black")
mtext("Non-Migrants",side=3,line=1,outer=F, font = 2)
mtext("% Agents",side=1,line=3,outer=F)
axis(2, at=num.sets, labels=plot.obj.names, las = 2)

dev.off()
